% 利用单平面光场进行Z轴配准

%  images path
DIR_path_LF_cell= ...
    { ...
    '/data8t/hryang/beads/20251103/SP/LFM/03_Reconstructed/'
    };

% list all images path
IMG_path_LF_cell={};
for i = 1:length(DIR_path_LF_cell)
    File_Path=DIR_path_LF_cell{i};
    File = dir(fullfile(File_Path,'*Z_1.tif'));
    IMG_path_LF_cell_i = cellfun(@(x) strcat(File_Path, x), {File.name}, 'UniformOutput',false);
    IMG_path_LF_cell_i = sort_nat(IMG_path_LF_cell_i);

    IMG_path_LF_cell=[IMG_path_LF_cell IMG_path_LF_cell_i];
end

N=length(IMG_path_LF_cell);
Results=zeros(N,1);

% process one by one
for i = 1:N

    LF_im_name = IMG_path_LF_cell{i};
    disp(['frame: ' num2str(i) '/' num2str(length(IMG_path_LF_cell)) '  ' LF_im_name]);

    [~, base, ~] = fileparts(LF_im_name);

    % ===== 读取 3D 体数据 (Y,X,Z) =====
    V = loadtiff(LF_im_name);
    % 统一转换为 double 计算（如需保持原类型保存，后面再 cast 回来）
    Vd = double(V);

    muPerSlice = squeeze(mean(mean(Vd, 1, 'omitnan'), 2, 'omitnan'));  % [Z×1]
    [muMax, zMax] = max(muPerSlice);
    Results(i)=zMax;

end

fprintf('完成：共处理 %d 个 3D TIFF，结果已保存到：%s\n', N, savePath);

%%
%对找到的beads的z轴坐标进行线性拟合
fit_xx=1:1:101;
fit_yy=flip(Results(:,1));

figure; plot(fit_yy);

fit_xx(34:end)=[];
fit_yy(34:end)=[];

p = polyfit(fit_xx,fit_yy,1);
z_fit = polyval(p, fit_xx);

figure;
plot(fit_xx, fit_yy,'o' ,fit_xx ,z_fit ,'-');
legend('data','linear fit');